Association of circulating gene expression signatures with stiffness following total knee arthroplasty for osteoarthritis: a pilot study

A subset of patients undergoing total knee arthroplasty (TKA) for knee osteoarthritis develop debilitating knee stiffness (reduced range of motion) for poorly understood reasons. Dysregulated inflammatory and immune responses to surgery correlate with reduced surgical outcomes, but the dysregulated gene signatures in patients with stiffness after TKA are poorly defined. As a consequence, we are limited in our ability to identify patients at risk of developing poor surgical outcomes and develop preventative approaches. In this pilot study we aimed to identify perioperative blood gene signatures in patients undergoing TKA for knee osteoarthritis and its association with early surgical outcomes, specifically knee range of motion. To do this, we integrated clinical outcomes collected at 6 weeks after surgery with transcriptomics analyses in blood samples collected immediately before surgery and at 24 h after surgery. We found that patients with stiffness at 6 weeks after surgery have a more variable and attenuated circulating gene expression response immediately after surgery. Our results suggest that patients with stiffness following TKA may have distinct gene expression signatures detectable in peripheral blood in the immediate postoperative period.

) isolated from PAXgene tubes collected at DOS and POD1 from 18 patients (8 cases and 10 matched controls) passed the RNA integrity quality controls and were used for RNA-seq, following standard procedures. A total of 100 ng of RNA were used to construct libraries, and sequencing was performed using an Illumina HiSeq 2500 at the WCM Epigenomics Core Facility using standard protocols 16 .
Bioinformatics analyses. After sequencing, the reads were processed using established pipelines at the David Z. Rosensweig Genomics Research Center at Hospital for Special Surgery (HSS). Briefly, reads were preprocessed using fastp 17 , which supports adapter trimming, low quality base trimming, and calculation of additional QC metrics. Trimmed, high quality reads were aligned to the target genome using STAR 18 . Low quality and multimapping alignments were filtered out using SAMtools 19 . Reads were counted within exons and summarized at the gene level using featureCounts 20 to produce a count matrix of reads-per-gene. These counts were analyzed for differential expression using the edgeR QLF framework 21 . Differentially expressed genes were used to perform QuSAGE pathway analyses 22 against MsigDB 23 pathways and gene sets. Transcription factor regulatory networks were generated from AnimalTFDB 3.0 24 and RegNetwork 25 , which include both transcription factors and transcriptional co-factors. All results were visualized with plotly 26 on an R Shiny 27 platform developed by the HSS Genomics Center.

Statistical analyses.
Analyses of demographics and range-of-motion at baseline and 6 weeks after TKA were conducted by the Department of Anesthesiology, Critical Care & Pain Management at HSS to identify patients with knee stiffness (cases) and matched non-stiff patients (controls). Samples were identified after enrollment end. Unpaired Student t-test was used to establish statistical significance between DOS and POD1 pathway and cell scores using GraphPad Prism 8 Software (GraphPad Software, San Diego, CA).

NanoString analyses identified surgery-induced gene expression signatures in PBMCs.
To establish circulating gene expression responses to TKA, first we analyzed RNA from PBMCs collected from 6 patients at DOS and POD1 using the human nCounter NanoString Immunology panel. This panel contains 594 genes grouped in different gene sets and permits the multiplexed evaluation of changes in gene expression associated with immune cell types and functional pathways. The initial assessment identified gene expression changes consistent with responses to surgery in all patients except one (Supplementary Figure S1). Chart review revealed that this patient received twice the dose of two anti-inflammatory medications (ketorolac and meloxicam) relative to the other patients. Thus, this patient's samples were excluded from the final analyses.   Fig. 2A. Differential gene expression analyses identified changes in S100A8, S100A9, SOCS3, CD163, CCR1, CCR2, CR1, BLC3, STAT3, and CXCL1 mRNA, as shown in the volcano plot representation in Fig. 2B. The top 20 upregulated genes (log2FC > 1.5, p < 0.05) on POD1 are shown in Table 1. We next used the differentially expressed genes (DEGs) identified by NanoString to perform pathway and cell score analyses. These analyses uncovered a relative enrichment in the immune response genes (including significant changes in XBP1, NCF4, GPI, BST1, or CEBPB mRNA; adjusted p value < 0.05), inflammatory response genes (with significant changes in S100A9, S100A8, CCR1, PTAFR, or CYBB; adjusted p value < 0.05), response to stress genes (with significant changes in S100A9, S100A8, PTAFR, MAPK14, or CYBB; adjusted p value < 0.05), and response to wounding genes (driven by increased expression of S100A9, S100A8 and CYBB; adjusted p value < 0.05) at POD1 (Fig. 2C). We also observed a relative enrichment in monocyte cell scores (driven by increased CD136 mRNA) accompanied by decreased T-cell scores (determined by changes in the expression of CD6, SH2D1A, CD3E and CD3D mRNA) at POD1 (Fig. 2D), consistent with reports showing increased monocytes and decreased T cells early after surgery 13 . See Supplementary Table S1 for the NanoString normalized counts and Supplementary  Table S2 for the NanoString pathway and cell-score analyses.  Table S3 for a summary of all DEGs. Comparison of the RNA-seq (PAXgene blood) and NanoString (PBMCs) datasets uncovered similar DEG profiles at POD1, including changes in CCR2, CD27, MARCO, PTGS2, RORC, S10018, S100A9, or SOCS3 (Fig. 3C). Together, these changes in gene expression indicate that the PAXgene RNA-seq data provide a valid screening tool to identify surgery-induced changes. In spite of the strong similarities between datasets, these comparisons also uncovered differences between the DEGs identified in PBMCs and the RNA-seq analyses in PAXgene blood RNA (Supplementary Figure S2), which can be driven by the different platforms used for analyses, the differential presence of polymorphonuclear leukocytes (PMNs), or the osmotic stress associated with Ficoll purification, but are also consistent with reports highlighting differences between RNA from purified PBMCs and PAXgene blood RNA 28,29 . Next, we performed QuSAGE pathway analyses 22 using RNA-seq data (Fig. 3D, FDR < 0.01). Comparison of POD1 relative to DOS samples uncovered changes in pathways associated with inflammatory and reparative responses, including increased expression of the IL1 (Fig. 3E), IFNG (Fig. 3F), and IL6 signaling (Fig. 3G), and decreased expression of the DNA repair pathway at POD1 (Fig. 3H). See Supplementary Table S4 for a summary of all differentially expressed pathways in POD1 vs. DOS samples. Taken together, our RNA-seq analyses in PAXgene samples confirmed the presence of surgery-induced signatures early after TKA, detectable in whole blood and characterized by an overall upregulation of inflammatory and stress responses consistent with our NanoString gene expression data in PBMCs.
Patients with knee stiffness at 6 weeks after TKA display distinct gene signatures at 24 h after surgery. Comparing cases and controls, we did not find significant differences for age (p = 0.95 by t-test), BMI (p = 0.32 by t-test), baseline flexion (p = 0.07 by t-test), baseline extension (p = 0.76 by t-test), baseline rangeof-motion (p = 0.09 by t-test), and sex distribution (p = 0.19 by Fisher's exact test). As expected, cases and controls had significant differences in post-operative range-of-motion (p = 0.0002 by t-test). Supplementary Table S5 summarizes the categorical demographics and the pre-and post-operative variables of cases and controls used for RNA-seq. All implants in the cases and controls included in this pilot study were cemented, with implant manufacturers similarly distributed between cases and controls. The same surgical approach was used for all patients (midline incision with a medial peripatellar arthrotomy). Further, we did not observe differences in tourniquet time between cases and controls (p = 0.94 by t-test).
We further analyzed our PAXgene RNA-seq dataset to identify perioperative gene signatures associated with stiffness at 6 weeks following TKA. We did not detect significant differences in gene expression comparing cases Table 1. Top 20 differentially expressed genes identified by NanoString analyses in PBMCs obtained the day-of-surgery (DOS) and at 24 h after surgery (POD1) from 5 patients undergoing total knee arthroplasty for knee osteoarthritis, with Log2 fold-change > 1.5 and p < 0.05. Log2 FC = Log2 fold change (POD1-vs-DOS). Lower = Lower confidence limit (log2). Upper = Upper confidence limit (log2).    4A) and 162 DEGs in cases (Fig. 4B). The Venn diagram representation in Fig. 4C depicts the unique and overlapping genes in cases and controls. We identified 111 overlapping DEGs for cases and controls, including ARG1, CCR2, CCR3, RETN, S100A9, S100A12, and TLR5 ( Fig. 4A-B), which represent a general response to TKA. Notably, the control group displayed 120 unique DEGs, including inflammatory mediators and damage-induced genes like IL1B, OSM, PTGS2, and S100A8 (Fig. 4D). The cases displayed 51 unique DEGs at POD1 relative to DOS, including CD36, GBP4, LYZ and MARCH1 (Fig. 4E). Thus, comparison of cases versus controls uncovered attenuated and more variable responses to surgery in cases, with 231 DEGs at POD1 in controls versus the 162 DEGs identified in cases. See Supplementary Table S6 for a summary of unique and overlapping genes in cases and controls. QuSAGE pathway analyses 22 using RNA-seq data also uncovered pronounced differences in cases and controls. Figure 5A summarizes the 97 differentially expressed signaling pathways identified in controls, comparing POD1 vs. DOS samples (FDR < 0.01), highlighting changes in the expression of the ATF2 pathway (Fig. 5B), pathways relevant to cytokines and inflammatory responses (Fig. 5C), the canonical NF-κB pathway (Fig. 5D), and the activated TLR4 signaling (Fig. 5E) after surgery. Similar to the network analyses using the combined POD1 and DOS samples (not shown), network analyses using the DEGs in the control group identified a transcription factor regulatory network involving GATA2, GADD45A, and CXCR3 (Fig. 5F). This transcriptional network was absent in the cases. Functional analyses using RNA-seq data from cases revealed a quite different transcriptional profile. QuSAGE analyses only uncovered 14 differentially expressed pathways in cases, comparing POD1 vs. DOS (FDR < 0.01, Fig. 5G). The majority of these pathways were common between cases and controls. The only 2 pathways unique to cases were the matrisome associated (Fig. 5H) and the steroid hormone biosynthesis (Fig. 5I). Important pathways involved in inflammatory, stress, and healing responses, including canonical NF-κB, p38  Taken together, our RNA-seq analyses confirmed the presence of surgery-induced circulating gene signatures following TKA and showed that changes in these early gene signatures are associated with stiffness at 6 weeks after surgery.

Discussion
In this pilot study, using comprehensive transcriptomics analyses in PBMCs and PAXgene whole blood RNA, we described the early circulating responses to surgery in patients undergoing TKA for OA. Integrating these analyses with clinical outcomes from our well characterized, prospectively enrolled cohort, we showed that changes in surgery-induced whole blood gene expression signatures associated with stiffness at 6 weeks after TKA can be detected early after surgery. Comparison of cases (patients who develop stiffness) versus controls uncovered attenuated and more variable responses to surgery in cases, with 231 DEGs at POD1 in controls versus the 162 DEGs identified in cases.
General responses to surgical and accidental trauma are well described 30,31 . Three recent studies, including our own, used targeted approaches and association with surgical outcomes to show changes in cytokine profiles and gene expression following TKA 9,14,15 . In the current study, using transcriptomics analyses in purified PBCMs www.nature.com/scientificreports/ and whole blood, we showed that TKA surgery leads to pronounced transcriptional changes in peripheral blood within 24 h, consistent with the previously described responses to trauma 30-32 and the changes described in patients recovering from total hip arthroplasty 13 . To address the known potential artifacts that can be introduced during sample collection, handling, storage, and processing of whole blood preservation systems 33 , we compared RNA-seq from PAXgene blood RNA with NanoString data from PBMCs collected from 6 consecutive patients. Integrating NanoString and RNA-seq datasets, we found that monocyte/macrophage gene signatures are enriched at 24 h after surgery, with a concomitant increase in the expression of alarmins (e.g., S100A8), pattern recognition receptor signatures (e.g., TLR4), and inflammatory cytokines (e.g., IL1B), and with decreased T-cell gene signatures. Macrophages exist within a phenotypic spectrum ranging from the classical (so-called M1) pro-inflammatory to the alternative (M2) state of activation, associated with the resolution of inflammation and healing. The coordinated action of these cellular components is critical to mediate host responses to tissue damage 34 . Our data do not permit single-cell resolution, so we cannot identify the relative enrichment of specific monocyte/macrophage subsets or activation states; however, we can identify a robust increase in monocyte signatures following surgery, including markers of M2-like macrophages like ARG1, CD163, and MARCO 35,36 . Soon after injury, and in parallel to the recruitment of macrophages and increased inflammation, there is also a compensatory phase characterized by reduced circulating T cell subsets, in part driven by increased migration into the periphery 13,37 . Our NanoString and RNAseq datasets agree with these previous reports and also show reduced T-cell gene expression at 24 h after surgery. Thus, the rapid and coordinated changes in gene expression signatures that we observed at 24 h after TKA surgery are consistent with the well-described general responses to injury, which are essential to drive optimal reparative processes.
Knee stiffness after TKA is often associated with the development of a fibrotic reaction in the knee joint (arthrofibrosis), indicating that impaired reparative responses lead to pathological scar formation concomitant with decreased range of motion [38][39][40][41] . In our study, comparison of cases (patients with stiff knees at 6 weeks after TKA) and controls (no stiffness) at baseline did not uncover significant differences in gene expression prior to surgery. However, comparison of the response to surgery in cases versus controls revealed a more variable and attenuated response to surgery in patients with stiffness after TKA. Comparing transcriptional profiles at a singlegene expression level we identified a common set of genes differentially expressed in response to surgery in both cases and controls, as well as subsets of genes that are unique and differentially expressed in controls (GATA2, IL1B, OSM, S100A8, and PTGS2) or in cases (LYZ, MARCH1, CD36, and RP2) after surgery.
When we evaluated functional gene signatures, we found that signaling pathways associated with stress, inflammation, and wound healing (including IL1, ATF2, or canonical NF-κB) are differentially expressed immediately after surgery only in controls, consistent with the expected response to injury and the immediate attempts to repair tissue damage 32,42 . The cases, however, did not display significant changes in the expression of these relevant signaling pathways at 24 h after surgery. This observation suggests that knee stiffness at 6 weeks post-TKA may be associated with changes in surgery-induced gene expression signatures that contribute to normal healing responses. Indeed, canonical NF-κB activation is observed in response to injury 43 , and the orchestration of canonical NF-κB targets is required for optimal tissue healing and repair and resolution of inflammation 44 . Thus, increased inflammatory and NF-κB-related responses in our control group can be interpreted as optimal reparative responses to trauma. However, we cannot conclude that the lack of changes in these pathways is proximate cause of stiffness in our cases. We did not evaluate whether the circulating/systemic responses mirror the inflammatory changes in the joint microenvironment that have been associated with fibrosis and stiffness 38,40 , and if the duration of these responses (or lack thereof) impacts outcomes, as is the case with persistent NF-κB activation leading to fibrotic healing 45 . Thus, the reproducibility, potential functional implications, and mechanistic relevance of our data in this context merit further investigation in larger and more diverse patient populations. Albeit not significant (p = 0.09 by t-test) the cases selected for this pilot study showed a relative enrichment in pre-operative stiffness (75% cases and 40% controls). While we did not observe differences in gene expression between cases and controls at baseline, we cannot rule out that (albeit not significant) the relative enrichment in pre-operative stiffness contributes to the attenuated gene expression observed in the cases after surgery. It is possible that the stimuli that lead to pre-operative stiffness modified the responses to surgery, akin to mechanisms whereby an initial exposure modify subsequent immune responses 46,47 . This should be evaluated in future work, comparing local and circulating signatures before and after surgery in patients without and with stiffness pre-operatively.
While the responses to surgery and the differences between cases and controls are robust and were obtained in a well-characterized patient cohort, our pilot study is not without limitations. Our data are correlative and do not allow us to establish direct mechanistic connections between peripheral gene expression and the development of stiffness and fibrosis after TKA, which has been shown to be associated with local inflammatory responses 38,40 . Our small sample size and the relative homogeneity of our patient population are limitations, and our data need to be validated and reproduced in larger and more heterogeneous cohorts. All patients enrolled in this study followed a standardized rehabilitation prototocol and we did not observe significant differences in baseline range-of-motion between cases and controls (p = 0.09). However, we did find that cases had higher pain Numeric Rating Scale (NRS) scores with movement prior to surgery (p = 0.02) despite no differences in baseline gene expression profiles. Preoperative knee pain predicts development of persistent pain following TKA 15,[48][49][50][51] ; however, the potential relationships between pain and other potential confounding and predisposing factors to knee stiffness, and their interaction with surgery-induced circulating responses, should be further analyzed in future studies and larger patient cohorts. The use of early stiffness (6 weeks) to identify cases and controls limits our ability to extrapolate our findings to patients who develop persistent or refractory knee stiffness and fibrosis requiring revision surgeries. However, 6 weeks after TKA is a clinically relevant time-point for determining the need for manipulation under anesthesia (MUA) as a means of recovering motion at the knee 11 www.nature.com/scientificreports/ of the cases in this pilot study (5 out of 8 patients who developed stiffness at 6 weeks after TKA) had stiffness requiring MUA. Finally, although PAXgene blood RNA tubes are a well-accepted and useful method for RNA collection and preservation in a clinical setting, RNA isolation from these samples is subject to technical and processing artifacts. These limitations, along with the small sample size used for RNA-seq and NanoString analyses, could have affected our ability to detect differences in cytokine transcripts with low expression levels and to establish comparisons with the protein data that we previously reported 9 .
In this pilot study, the changes in gene expression associated with stiffness after TKA were obtained in samples collected within 24 h after surgery. Future studies should aim to integrate these perioperative signatures with changes in gene expression detected at later time-points to identity signatures associated with the development of refractory knee stiffness. Future work should also aim to comprehensively address the predisposing factors leading to the development of stiffness and fibrosis, the association of stiffness with abnormal local and circulating immune responses, and the patient-intrinsic factors that contribute to the variable responses to treatment. These future studies should establish clinically relevant correlations between gene expression signatures and patient outcomes, aiming to identify patients and risk and prevent suboptimal outcomes. Additional studies that integrate clinical, histological, and multimodal cellular and genomics analyses on tissues retrieved at the time of primary TKA and revision surgeries can provide information about the pathways that, when dysregulated, contribute to poor recovery from TKA and the development of knee stiffness after surgery.
In conclusion, our results show that peripheral gene signatures can be used to evaluate pathways involved in the responses to surgery and that patients with stiffness following TKA may have dysregulated gene signatures detectable in the acute postoperative period. Notably, these gene expression signatures were detected from whole blood samples collected in PAXgene blood RNA tubes for transcriptomic analyses, suggesting a clinically feasible approach to developing molecular methods for predicting patients at risk of developing complications following surgery.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request. The RNA-seq sequencing data have been deposited at the database of Genotypes and Phenotypes (dbGaP) under dbGaP accession code phs002927.v1.p1.